home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / tmth.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-14  |  6.7 KB  |  217 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5.  
  6. #include "newmatap.h"
  7. //#include "newmatio.h"
  8.  
  9. void Print(const Matrix& X);
  10. void Print(const UpperTriangularMatrix& X);
  11. void Print(const DiagonalMatrix& X);
  12. void Print(const SymmetricMatrix& X);
  13. void Print(const LowerTriangularMatrix& X);
  14.  
  15. void Clean(Matrix&, Real);
  16.  
  17.  
  18.  
  19.  
  20. void trymath()
  21. {
  22. //   cout << "\nSeventeenth test of Matrix package\n";
  23. //   cout << "\n";
  24.    Tracer et("Seventeenth test of Matrix package");
  25.    Exception::PrintTrace(TRUE);
  26.  
  27.  
  28.    {
  29.       Tracer et1("Stage 1");
  30.       int i, j;
  31.       BandMatrix B(8,3,1);
  32.       for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
  33.      { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
  34.  
  35.       DiagonalMatrix I(8); I = 1; BandMatrix B1; B1 = I;
  36.       UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
  37.       Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
  38.       Matrix A = B; BandMatrix C; C = B;
  39.       Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
  40.  
  41.       C.ReDimension(8,2,2); C = 0.0; C.Inject(B);
  42.       Matrix X = A.t();
  43.       B1.ReDimension(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
  44.  
  45.       Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
  46.       BandLUMatrix BLU = B.t();
  47.       BI = BLU.i(); A = X.i()-BI; Clean(A,0.000000001); Print(A);
  48.       X.ReDimension(1,1);
  49.       X(1,1) = BLU.LogDeterminant().Value()-B.LogDeterminant().Value();
  50.       Clean(X,0.000000001); Print(X);
  51.  
  52.       UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
  53.       DiagonalMatrix D; D << B;
  54.       Print( Matrix(L + (U - D - B)) );
  55.  
  56.  
  57.  
  58.       for (i=1; i<=8; i++)  A.Column(i) << B.Column(i);
  59.       Print(Matrix(A-B));
  60.    }
  61.    {
  62.       Tracer et1("Stage 2");
  63.       BandMatrix A(7,2,2);
  64.       int i,j;
  65.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  66.       {
  67.      int k=i-j; if (k<0) k = -k;
  68.      if (k==0) A(i,j)=6;
  69.      else if (k==1) A(i,j) = -4;
  70.      else if (k==2) A(i,j) = 1;
  71.      A(1,1) = A(7,7) = 5;
  72.       }
  73.       DiagonalMatrix D(7); D = 0.0; A = A - D;
  74.       BandLUMatrix B(A); Matrix M = A;
  75.       ColumnVector V(3);
  76.       V(1) = LogDeterminant(B).Value();
  77.       V(2) = LogDeterminant(A).Value();
  78.       V(3) = LogDeterminant(M).Value();
  79.       V = V / 64 - 1; Clean(V,0.000000001); Print(V);
  80.       ColumnVector X(7);
  81.  
  82. #ifdef ATandT
  83.       Real a[7];
  84.       // the previous statement causes a core dump in tmti.cpp
  85.       // on the HP9000 - seems very strange. Possibly the exception
  86.       // mechanism is failing to track the stack correctly. If you get
  87.       // this problem replace by the following statement.
  88. //      Real* a = new Real [7]; if (!a) exit(1);
  89.       a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  90. #else
  91.       Real a[] = {1,2,3,4,5,6,7};
  92. #endif
  93.       X << a;
  94. #ifdef ATandT
  95.       delete [] a;
  96. #endif
  97.       M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
  98.       Clean(M,0.000000001); Print(M);
  99.  
  100.  
  101.       BandMatrix P(80,2,5); ColumnVector CX(80);
  102.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  103.       { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
  104.       for (i=1; i<=80; i++)  CX(i) = i*100.0;
  105.       Matrix MP = P;
  106.       ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
  107.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  108.  
  109.       V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
  110.       RowVector XX(1);
  111.       XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
  112.       Clean(XX,0.000000001); Print(XX);
  113.  
  114.       LowerBandMatrix LP(80,5);
  115.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  116.       { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
  117.       MP = LP;
  118.       XX.ReDimension(4);
  119.       XX(1) = LogDeterminant(LP).Value();
  120.       XX(2) = LogDeterminant(MP).Value();
  121.       V1 = LP.i() * CX; V2 = MP.i() * CX;
  122.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  123.  
  124.       UpperBandMatrix UP(80,4);
  125.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  126.       { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
  127.       MP = UP;
  128.       XX(3) = LogDeterminant(UP).Value();
  129.       XX(4) = LogDeterminant(MP).Value();
  130.       V1 = UP.i() * CX; V2 = MP.i() * CX;
  131.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  132.       XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
  133.    }
  134.    {
  135.       Tracer et1("Stage 3");
  136.       SymmetricBandMatrix SA(8,5);
  137.       int i,j;
  138.       for (i=1; i<=8; i++) for (j=1; j<=8; j++)
  139.      if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
  140.       DiagonalMatrix D(8); D = 10; SA = SA + D;
  141.  
  142.       Matrix MA1(8,8); Matrix MA2(8,8);
  143.       for (i=1; i<=8; i++)
  144.      { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
  145.       Print(Matrix(MA1-MA2));
  146.  
  147.       D = 10; SA = SA.t() + D; MA1 = MA1 + D;
  148.       Print(Matrix(MA1-SA));
  149.  
  150.       UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
  151.       D << SA; UB << SA; LB << SA;
  152.       SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
  153.       BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
  154.  
  155.       SymmetricBandMatrix A(7,2); A = 100.0;
  156.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  157.       {
  158.      int k=i-j;
  159.      if (k==0) A(i,j)=6;
  160.      else if (k==1) A(i,j) = -4;
  161.      else if (k==2) A(i,j) = 1;
  162.      A(1,1) = A(7,7) = 5;
  163.       }
  164.       BandLUMatrix C(A); Matrix M = A;
  165.       ColumnVector X(8);
  166.       X(1) = LogDeterminant(C).Value() - 64;
  167.       X(2) = LogDeterminant(A).Value() - 64;
  168.       X(3) = LogDeterminant(M).Value() - 64;
  169.       X(4) = SumSquare(M) - SumSquare(A);
  170.       X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
  171.       X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
  172.       X(7) = Trace(M) - Trace(A);
  173.       X(8) = Sum(M) - Sum(A);
  174.       Clean(X,0.000000001); Print(X);
  175.  
  176. #ifdef ATandT
  177.       Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  178. #else
  179.       Real a[] = {1,2,3,4,5,6,7};
  180. #endif
  181.       X.ReDimension(7);
  182.       X << a;
  183.       X = M.i()*X - C.i()*X * 2 + A.i()*X;
  184.       Clean(X,0.000000001); Print(X);
  185.  
  186.  
  187.       LB << A; UB << A; D << A;
  188.       BandMatrix XA = LB + (UB - D);
  189.       Print(Matrix(XA - A));
  190.  
  191.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  192.       {
  193.      int k=i-j;
  194.      if (k==0) A(i,j)=6;
  195.      else if (k==1) A(i,j) = -4;
  196.      else if (k==2) A(i,j) = 1;
  197.      A(1,1) = A(7,7) = 5;
  198.       }
  199.       D = 1;
  200.  
  201.       M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
  202.       M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
  203.       M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
  204. #ifdef __GNUG__
  205.       Matrix MUB = UB; Matrix MLB = LB;
  206.       M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
  207.       M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
  208. #else
  209.       M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
  210.       M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
  211. #endif
  212.    }
  213.  
  214.  
  215. //   cout << "\nEnd of Seventeenth test\n";
  216. }
  217.